## ----------------------------------------
# This file contains R functions that can be re-used for many analyses
# It loads common libraries that are frequently used.
# Written on: 03/21/2014
# Last update: 09/06/2014
# Author: Thong Nguyen
## ----------------------------------------

options(stringsAsFactors=FALSE)
options(check.names=FALSE)
suppressMessages(library(parallel))
suppressMessages(library(plyr))
suppressMessages(library(data.table))
suppressMessages(library(RMySQL))
suppressMessages(library(sqldf))
suppressMessages(library(plotrix))
suppressMessages(library(bear))
suppressMessages(library(ggplot2))
suppressMessages(library(rmarkdown))

my.ucsc_connect <- function() {
	dbConnect(MySQL(), user = "genome", dbname = "hg19", host = "genome-mysql.cse.ucsc.edu")
}


my.aplle_connect <- function() {
	dbConnect(MySQL(), user = "aplle_owner", password="ngsequencing", dbname = "aplle", host = "mysqlprd.gene.com", port = 3320)
}

my.resthongn_connect <- function() {
	dbConnect(MySQL(), user = "resthongn_owner", password="th139n", dbname = "resthongn", host = "mysqld01.gene.com", port = 5606)
}
my.thong_connect <- function(host = "thongn.wifi.gene.com") {
	dbConnect(MySQL(), user = "thong", password="piraat103", dbname = "gene", host = host)
}
my.load_vcf <- function(var_file) {
	check_file <- function(x) {
		l = readLines(x)
		l = substr(l, 1, 2)
		l1 = l[l=="##"]
		c(length(l1), length(l) - length(l1))
	}
    cf = check_file(var_file)
    if (cf[2] < 2) return (NA)

    outdt = read.table(var_file, sep = "\t", header = TRUE, skip = cf[1], comment.char = "")
	colnames(outdt)[1] = "varid"
	outdt
}
mnlapply <- function(X, FUN, cores = 4, queue = "veryshort", vars = NULL) {

	# Create R run script
	if (is.null(vars)) {
		save.image("X.rda", compress = FALSE)
	} else {
		save(file = "X.rda", list = vars, compress = FALSE)
	}

	tmp = "./tmpdata"
	dir.create(tmp)
	
	run_script = sprintf("%s/run.R", tmp)
	system(sprintf("Rheader %s", run_script))
	
	sub_script= "	
				load(\"X.rda\")
				cmd_args = commandArgs(TRUE)
				FUN = cmd_args[1]
				x = cmd_args[2]
				rdsfile = cmd_args[3]
				print(x)
				out = do.call(FUN, list(x))
				saveRDS(out, file = rdsfile, compress = FALSE)
				"
	write(sub_script, file = run_script, append = TRUE) 
	
	# Run script for each item in X
	execute <- function(i) {
		x = X[[i]]
		job = sprintf("%s/job%s", tmp, i)
		rdsfile = paste(job, ".rds", sep = "")
		
		command = sprintf("%s %s %s %s", run_script, FUN, x, rdsfile)
		bsub = sprintf( "bsub -q %s -sp 25 -n %s -K -R \"span[hosts=1]\" -e %s.err -o %s.out \"%s\"", 
					    queue, cores, job, job, command )
		system(bsub)
		readRDS(rdsfile)
	}
	
	njobs = length(X)
	mclapply(1:njobs, execute, mc.cores = njobs)
}


## Dump heteroplasmy table from MySQl DB to VCF format
my.vcf_dump <- function() {
source("http://salthill.googlecode.com/svn/trunk/Rlib.R")

con <- dbConnect(MySQL(), user = "thong", password="piraat103", dbname = "gene", host = "thongn.wifi.gene.com")

tabnames = c("392_heteroplasmy", "115_heteroplasmy", "488_heteroplasmy2", "196_heteroplasmy")
bam_paths = c("/gnet/is3/research/data/dnaseq/analysis/aplle/technology/rnaseq/392_NCCRCC_mito/WholeMito/clipOverlap",
			  "/gnet/is3/research/data/dnaseq/analysis/aplle/technology/rnaseq/115_Colon/WholeMito/clipOverlap",
			  "/gnet/is3/research/data/dnaseq/analysis/aplle/technology/whole_genome/488_NCCRCC_mito/WholeMitoSE/BAMs",
			  "/gnet/is3/research/data/dnaseq/analysis/aplle/technology/exome/196_colon_mito/WholeMito/clipOverlap")
names(bam_paths) = tabnames

create_vcf <- function(tabname, samtype) {
	bam_path = bam_paths[tabname]
	dir.create(tabname)
	sql = sprintf("select Chrom as `#CHROM`, 
				  Pos as POS, 
				  varid as ID, 
				  Ref as REF, 
				  Alt as ALT,
				  1000 as QUAL, 
				  'PASS' as FILTER, 
				  total%s as INFO, 
				  AltFreq%s as FORMAT, 
				  SamID%s as SAMPLE 
				  from %s where AltFreq%s > 0.05 and AltFreq%s < 0.95 and total%s > 100", 
				  samtype, samtype, samtype, tabname, samtype, samtype, samtype)

	tab = dbGetQuery(con, sql)
	sams = unique(tab$SAMPLE)

	write_vcf <- function (sam) {
		x = subset(tab, SAMPLE == sam)
		bam = sprintf("%s/%s_16_1.bam", bam_path, sam)
		vcf = sprintf("/gnet/is3/research/data/dnaseq/analysis/thongn/projects/Mito_488/analysis/heteroplasmy2/vcf/%s/%s.vcf", tabname, sam)
		write.table(x = x, file = vcf, sep = "\t", quote = F, row.names = F, col.names = T)
		data.frame(SamID = sam, BAM = bam, vcf =  vcf)
	}
	data.frame( rbindlist( mclapply(sams, write_vcf, mc.cores = length(sams)) ) )
}

tabnames = c("392_heteroplasmy", "115_heteroplasmy", "488_heteroplasmy2", "196_heteroplasmy")

df1 = data.frame( rbindlist( lapply(tabnames, create_vcf, samtype = "") ) )
df2 = data.frame( rbindlist( lapply(tabnames, create_vcf, samtype = "_normal") ) )

df = rbind(df1, df2)

write.table(x = df, file = "vcf.txt", sep = "\t", quote = F, col.names=T, row.names=F)

dbDisconnect(con)
}


## New function here
